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Abstract 

The Levy, jumping process, defined in terms of the jumping size distribution 
and the waiting time distribution, is considered. The jumping rate depends on 
the process value. The fractional diffusion equation, which contains the variable 
diffusion coefficient, is solved in the diffusion limit. That solution resolves itself 
to the stretched Gaussian when the order parameter fi ^ 2. The truncation of 
the Levy flights, in the exponential and power-law form, is introduced and the 
corresponding random walk process is simulated by the Monte Carlo method. The 
stretched Gaussian tails are found in both cases. The time which is needed to reach 
the limiting distribution strongly depends on the jumping rate parameter. When 
the cutoff function falls slowly, the tail of the distribution appears to be algebraic. 
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1 Introduction 



Transport in pliysical systems can be described in terms of a jumping process 
which is completely defined by two probabihty distributions. They determine 
spatial and temporal characteristics of the system and are called the jumping 
size, X{x), and waiting time, w{t), probability distributions, respectively. The 
stochastic trajectory consists of infinitely fast jumps which are separated by 
intervals when the Brownian particle is at rest (e.g., due to traps). Usually, 
the distributions A and w are regarded as mutually independent and processes 
which they generate are studied in the framework of the decoupled continuous 
time random walk (CTRW). If the distribution w{t) is Poissonian, the jumps 
take place at random, uniformly distributed times and the corresponding pro- 
cess is Markovian. The algebraic form of w{t), which possesses long tails, is 
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of particular interest; it leads to long rests and then to a weak, subdiffusive 
transport [1]. The jumping size distribution A can assume - in order to be sta- 
ble - either the Gaussian form, or obey the general Levy distribution which 
corresponds to processes with long tails, frequently encountered in many areas 
of physics, sociology, medicine [2] and many others. In the diffusion limit of 
small wave numbers, the master equation for the decoupled CTRW, and for 
the case of the Levy distributed jumping size, resolves itself to the fractional 
diffusion equation with a constant diffusion coefficient. 

However, many physical phenomena require taking into account that the space 
has some structure and that the diffusion coefficient must actually be variable. 
It is the case when one considers the transport in porous, inhomogeneous me- 
dia and in plasmas, as well as for the Levy flights in systems with an external 
potential [3] . In the Hamiltonian systems, the speed of transport depends on 
regular structures in the phase space [4] . Any description of the diffusion on 
fractals must involve the variable diffusion coefficient, and that one in the 
power-law form [5,6]. Similarly, the diffusion on multifractal structures can be 
regarded as a superposition of the solutions which correspond to the individ- 
ual fractals [7]. The power-law form of the diffusion coefficient has been also 
used to describe e.g. the transport of fast electrons in a hot plasma [8] and 
the turbulent two-particle diffusion [9]. The presence of long jumps indicates a 
high complexity and the existence of long-range correlations. One can expect 
that especially in such systems the waiting time depends on the position. The 
jumping process, stationary and Markovian, which takes into account that 
dependence, has been proposed in Ref.[10]: the distribution w{t) is Poissonian 
with a a;-dependent jumping rate In this paper, we consider that process 
for the Levy distributed X{x) and solve the corresponding fractional equa- 
tion. The solution consistently takes into account that the symmetric Levy 
distribution, defined by the characteristic function exp(— < < 2, 
exhibits two qualitatively different kinds of behaviour in its asymptotic limit: 
it is either algebraic, ~ for < < 2 or Gaussian for fi — 2. 

The Levy process is characterized by very long jumps since the tails of the 
Levy distribution fall slowly and the second moment is divergent. In a concrete, 
realistic problem, however, the available space is finite and the asymptotics 
differs from the Levy tail. One can take that into account by introducing a 
cutoff to the stochastic equations by multiplying the jumping size distribution 
- in the Levy form - by a function which falls not slower than usually 
by the exponential or the algebraic function. Since the resulting probability 
distribution has the finite variance, in the diffusion approximation the process 
is equivalent to the Gaussian limit /i ^ 2. 

Restrictions on the jump size are necessary also for problems connected with 
the random transport of massive particles. Since the velocity of propagation 
is then finite, the jump, which takes place within a given time interval, must 
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have a finite length. If the jump size is governed by the particle velocity, a 
coupling between spatial and temporal characteristics of the system emerges 
in the framework of CTRW. That coupled form of the CTRW is known as 
the Levy walk [11,12,13]. It has many applications, e.g. for the turbulence [14] 
and the chaotic diffusion in Josephson junctions [11]. 

The aim of this paper is to study the limit in which the algebraic tail of the 
Levy process becomes the exponential function. Utilizing the result for — > 2 
can improve the solution of the fractional equation itself and it is suited to 
describe the truncated Levy flights in the diffusion approximation. The paper 
is organized as follows. In Sec. II we solve the fractional equation and derive the 
stretched Gaussian asymptotics in the limit fi ^ 2. Sec. Ill is devoted to the 
truncated Levy flights: we present the probability distributions, obtained from 
simulations of random walk trajectories, for both exponential and power-law 
form of the cutoff. The results are summarized in Sec.IV. 



2 The stretched-Gaussian limit of the Levy process 

We consider the random walk process defined by the waiting time probability 
distribution w{t) and the jump-size distribution X{x). They are of the form 

w{t) = i/(x)e-^^^)*, (1) 
where — \x\~^ {9 > —1) and 

oo 

X{x) = ^/2/^ J ex.pi-K^'k") cos{kx)dk, (2) 





respectively. The latter expression corresponds to the symmetric Levy distri- 
bution. The master equation for the above process is of the form [10] 

d f 

—p(x,t) = —iy(x)p(x,t)+ / iy(x')\(\x — x'\)p(x' ,t)dx' . (3) 
at J 



In the diffusion limit of small wave numbers, the Eq.(3) can be reduced to the 
fractional equation. Indeed, the Fourier transform of jump-size distribution, 
X{k) — exp(— can be expanded \{k) = 1 — K^\k\^^ + ... and the 
master equation in the Fourier space becomes the following equation [15] 

^ -K^\kmu(x)p(x,t)]. (4) 
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The above approximation means that the summation over jumps is substi- 
tuted by the integral and it agrees with the exact result if the jumps are 
small, compared to the entire trajectory. One can demonstrate by estimating 
the neglected terms in the Euler-Mclaurin summation formula that the ap- 
proximation fails near the origin and for < 1 [16]. In the other cases, the 
solution of the fractional equation converges with time to the exact result. 
Therefore, in the following, we assume /i > 1. 

The inversion of the Fourier transforms in Eq. (4) yields the fractional equation 



where d^/d\x\^ is the Riesz-Weyl derivative, defined for 1 < < 2 in the 
following way [17] 



^\x\^' 2cos(7r/i/2)r(2 — /i) (9x2 j — x' 



-oo 



Eq. (5) for the constant diffusion coefficient - often generahzed to take into ac- 
count non-Markovian features of the trapping mechanism in the framework of 
CTRW by substituting the simple time derivative by the fractional Riemann- 
Liouville derivative [1] - is frequently considered and solved by means of a 
variety of methods. The solution resolves itself to the Levy stable distribution 
with the asymptotic power-law x-dependence and divergent second moment. 
The method of solution which is especially interesting for our considerations 
involves the Fox functions. The well-known result of Schneider [18] states that 
any Levy distribution, both symmetric and asymmetric, can be expressed as 
In order to solve the Eq.(5) for the case of the variable jumping rate, 
= we conjecture that the solution also belongs to the class of func- 
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v{x) 

tions hI'I and it has the scahng property. Then the ansatz is the following 



p{x,t) = Na{t)Hl^ 



a{t)\x\ 



(ai, Ai), (02,^2) 



(61, Bi), (62,^2) 



(7) 



where N is the normalization constant. The Fox function is defined in the 
following way [19,20] 
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where 
X(s) 



nrrfe-i?,s)mr(i-a, + A,g) 
r(i - + B,s) Yt+i r(a, - A,s) • 



(9) 



The contour L separates the poles belonging to the two groups of the gamma 
function in the Eq.(9). Evaluation of the residues leads to the well-known 
series expansion of the Fox function: 
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We will try to solve the fractional equation (5) by inserting the function (7). 
This procedure, if successful, would allow us to find conditions for the coeffi- 
cients and the function a{t). The idea to assume the solution in the form (7) 
is motivated by the following property of the Fox function 
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which eliminates the algebraic factor by means of a simple shift of the coef- 
ficients. One can demonstrate that Eq.(5) cannot be satisfied, in general, by 
the function (7) for any choice of the parameters. However, as long as we are 
interested only in the diffusion limit of small wave numbers (large the 
higher terms in the characteristic function expansion can be neglected. There- 
fore we require that the Eq.(5) should be satisfied by a function which agrees 
with the exact solution only up the terms of the order | in the Fourier space. 
Note that this condition does not introduce any additional idealization since 
the Eq.(5) itself has been constructed on the same assumption: the higher 
terms in the |A;| expansion of X{x) have been also neglected. 

First, we need the Fourier transform of the Fox function. Since the process is 
symmetric, we can utilize the formula for the cosine transform which yields 
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also a Fox function but of the enhanced order: 
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The function p{k,t) is then of the order H^'l. The Fourier transform of the 
function pe — x~^p{x, t) can be obtained in the same way. Next, we insert the 
appropriate functions into the Eq.(4) and expand both sides of the equation by 
using the formula (10). Let us denote the expansion coefficients of the functions 
p{k,t) and peik.t) by h^^^, and /i^^^, respectively; a assumes the values 1 and 
2. Applying the Eq.(lO) yields 

p{k,t) = /ii,o + + h2,oH^~'''^'^'-^ + h2,i\kf~"'^'^"~^ + K^k^ + • • • (13) 



and 

Pe{K t) = h^l + h^^l\k\ + h^l\k\^'-^^-^'^^y^^~' + . . . . (14) 



After inserting the above expressions to the Eq. (4) , we find some simple rela- 
tions among the coefficients of the Fox function by comparison of the expo- 
nents. To get the term |A;|^ on Ihs, which corresponds to the term k^ on rhs, we 
need the condition (2 — ai)/Ai — 1 = /x. Moreover, we attach the third term on 
rhs to the first one by putting {l — ai-\-dAi)/Ai — l — 0; it is not possible to bal- 
ance the third term by another one on the Ihs. The above conditions determine 
two coefficients of the Fox function: Oi = 1 — — 9) / {ix+9) and >li = l/{ix+9). 

(0) 

The coefficients hi i and h\ [ vanish identically since the gamma function in 
the denominators has its argument equal 0. The only remaining term can be 
eliminated by assuming the condition 1 — 62 — -62(1 — ^) = 0; then /i2,o = 
since that term contains the function r(l — &2 — -62(1 — ^)) in the denominator. 
Therefore, such choice of the coefficients reduces the function (13) to the two 
terms and it makes it identical with the probability distribution for the Levy 
process in the diffusion limit. Finally, the Eq.(4) becomes a simple differential 
equation which determines the function a{t)\ 

i = K^^r'^^, (15) 

'12,1 
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where h^^ — h^iQ + /i2^o ^^<i C(0 = ^ The solution 

-i/in+e) 

a{t)^ 




(16) 



corresponds to the initial condition p{x, 0) = 5{x). The coefficient /i2,i can be 
determined directly from Eq.(lO), whereas h^^ — (27r)~^a~^p(0, t) — 7r~^a~^ Io° p(x, t)dx 
can be expressed in terms of the Mellin transform from the Fox function, x(s), 
and then easily evaluated. 

In the limit — > 2, the asymptotic behaviour of the fractional equation 
changes qualitatively. It is no longer algebraic; the tails of the distribution, 

and consequently the tails of the Fox function, have to assume the exponen- 
tial form for = 2. It is possible only if the algebraic contribution from the 
residues vanishes. If n = 0, all poles are outside the contour L and H^^^ is 
given by the integral over a vertical straight line: 



w+ioo 

^r.° = ^ / X{s)z^ds, (17) 

w—ioo 

where w < Re{bi/Bi). The above integral is usually neglected if n > because 
it is small compared to the contribution from the residues. 

The Fox function in the required form can be obtained from (7) by applying 
the reduction formula if the coefficients in the main diagonal are equal. This 
demand imposes an additional condition on 62 and B2 which allows us to 
determine these coefficients: 62 = 1 - (1 - 6')/(2 + 6') and B2 = 1/(2 + 6). 
Note that the above choice of (62,-82) yields hi^2 = in Eq.(13) and then the 
solution of Eq.(5) is exact up to the order k"^ for any /x. The most general 
solution of the fractional equation in the form of the function H22, which 
involves the required conditions, is the following 



p{x,t) ^ Na{t)Hl'^ 



a{t)\x\ 



1 ~ fel' («2,^) 



(61, Bi), (1 - 2^) 



(18) 



The coefficients in the main diagonal have a simple interpretation. Applying 
Eq.(lO) to (18) reveals that p{x,t) behaves as Ixj^^/'^^ for |x| and then 
the parameters bi and Bi are responsible for the shape of the probability 
distribution near the origin. On the other hand, the parameters Oi and Ai 
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determine the asymptotic shape of the distribution. It can be demonstrated 
by applying the following property of the Fox function 
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and by expansion according to Eq.(lO): the leading term is of the formp(a;, t) 



\x 



(2-ai)/Ai 



\X\ 



oo). Now it becomes clear why our method 
of solving the Eq.(5) did not determine the parameters {bi,Bi). Since we 
neglected higher terms in the A;-expansion, the region of small |a;| remained 
beyond the scope of the approximation. However, we can supplement that so- 
lution by referring directly to the master equation (3) which reveals a simple 
behaviour near the origin. That equation is satisfied by the stationary solu- 
tion l/iy{x) — \x\^, for any normalizable X{x). Obviously, such p{x,t) cannot 
be interpreted as the probability density distribution since the normalization 
integral diverges in infinity but it properly reproduces that distribution for 
small \x\. Therefore we obtain the additional condition 6i = 6Bi which im- 
proves the agreement of our solution with the solution of the master equation. 
The probabihty distribution for the case ^ = can be easily found by the di- 
rect solution of the fractional equation which is exact and yields the following 
values of the parameters: bi — 0, Bi — 1, a2 — 1/2, and A2 — 1/2. 

In the case = 2, the main diagonal in Eq.(18) can be eliminated and the 
solution of the fractional equation (5) in the limit — > 2 takes the form 



p(x,t) ^ Na(t)Hl'" 





{a2, A2) 


a{t)\x\ 






{bi,B,) _ 



(20) 



The function a{t) is given by Eq.(16) in the following form 

-i/{2+e) 



a{t) 



K\2 + e)^t 
h2 



(21) 



where Hq = r(6i + Bi(l-^))/r(a2+^2(l-^)) and h2 = T{bi+3Bi)/T{a2+3A2) 
are the coefficients of the /c-expansion of the functions \x\^^p{x,t) and p{x,t), 
respectively. The asymptotic expression for p{x, t) can be obtained from the 
estimation of the integral (17) by the method of steepest descents [21]. The 
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result reads 

Hli{z) ^ e*"("'-^/2)^(^e'"'^), (22) 



where 



-1 oo 

E{z) = -— 53Cj(/3a"^)(^-"'-^')/"exp(/3a"^)^/" (23) 



and a ^ Bi- A^, a' = - 6i + 1/2, = A^'jB^\ z = a\x\. The final 
result appears to be a stretched Gaussian, modified by an algebraic factor and 
a series which converges to a constant for 2; — > oo: 

1 00 
Hl'^{z) ^ exp(-/3^/"a^^/")/3(^-"')/"Q;^-"' J2yiyCjp-^/"a-h-- 



The above expression has been obtained by Wyss [22] as the expansion of 
-^2,3 (^) which satisfies a generalized diffusion equation. It is the integral equa- 
tion in respect to the time variable (non-Markovian) which resolves itself to 
the fractional diffusion equation; that equation is commonly used to handle 
the subdiffusive processes in the framework of the CTRW [12,23.1]. The func- 
tion which contains the stretched Gaussian, modified by the power-law factor, 
is used as the asymptotic form of the propagator for diffusion on fractals, e.g. 
on the Sierpihski gasket [24]. 

The coefficients Cj are defined by means of the following expression 
r(l-a. + A..) ^ g c, 



They can be explicitly evaluated by a subtraction of the consecutive terms 
and by taking the limit s — > 00. More precisely, Cj are given by the following 
recurrence formula 



d = lim 



GT{as + a')-Co- 



Ci Cj 



as + a' {as + a') . . . {as + a' + j — 1) 

X {as + a'){as + a' + 1) . . . {as + a' + j - 1) (2 

where 

Co = v^a-'-'/^Al^'-^'Bl'-'/^ (27) 
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has been obtained from the expansion of the gamma functions by means of the 
Stirhng formula. The exponent of the stretched Gaussian, is connected 
with higher moments of the distribution p{x, t) and it cannot be determined 
in the framework of the diffusion approximation. 

All moments of the distribution p{x,t) are convergent. In particular, the vari- 
ance, which determines the diffusion properties of the system, is given by the 
expansion coefficient /i2 in a simple way: 

(^') = -^P(^^ = ^20-' ~ . (28) 



For 9 = the diffusion coefficient D = limt^^{x'^) (t) /2t assumes a finite value 
and the diffusion is normal. The case 9^0 means the anomalous diffusion: 
either the enhanced one for ^ < 0, or the subdiffusion for ^ > 0; the diffusion 
coefficients are then oo or 0, respectively. The kind of diffusion depends only 
on 9 and it is not sensitive on free parameters. The anomalous diffusion is 
frequently encountered in physical phenomena, in particular in complex and 
disordered systems [25], as well as in dynamical systems [4]. 

On the other hand, the fractional equation (5) for /x = 2 assumes a form of 
the simple diffusion equation 

-^ = ^ d^^ ^'^^ 



which can be solved exactly just by assuming the scaling form of the solution 
p{x,t) = a{t)f{a{t)x). The functions a{t) and f{ax) are derived by inserting 
to the Eq.(29) and by separation of the variables [26]. That procedure finally 
yields 

p{x,t) = N{KH)-m\x\'e^p {- J{I+9)h ) " 



Eq.(29) follows from the master equation (3) with the Gaussian X{x), when one 
neglects all terms higher than the second one in the Kramers- Moyal expansion. 
That procedure is justified if jumps are small and ^{x) is a smooth function 
[27]. We can expect that the solution of the master equation converges with 
time to Eq.(30) and this convergence is fast for small \9\. Gonvergence of the 
tails must be slow, because the contribution from large jumps is substantial 
for large jxj. Indeed, we will demonstrate in the following that for large |^|, in 
particular for 9 close to —1, a very long time is required. 

The result (30) is useful for further improvement of the solution (18). The 
comparison of Eqs. (24) and (30) yields the conditions for the Fox function 
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coefficients which ensure the proper hmit /i 2: hm^^2(-Bi —^2) = 1/(2 + ^) 
and hm^_>2(l — — ^- By inserting the coefficients which follow from those 
limiting values to the Eq.(18) and by assuming, in addition, that Bi = 1, we 
obtain a particular solution of the fractional equation (5) in the form 



p{x,t) ^ Na{t)Hlll 



a{t)\x\ 



i-e 1 \ 



+ 



e{i+e) 
2+e 



,1 



2+e) 



{^1 1)) (1 2+61 2+e' 



■ (31) 



Therefore, the above solution of the fractional equation takes into account the 
Kramers-Moyal result (30) in the limit — > 2 and its behaviour in the origin 
agrees with the master equation. The hmit ^ — > corresponds to the Levy 
process. 



3 Truncated Levy flights 



The Levy flights are characterized by the power-law tail with the exponent 
smaller than 3 which implies the infinite variance. However, in physical sys- 
tems, which are limited in space, the variance must always be finite. The 
finiteness of the available space must be taken into account in any attempt 
to simulate the random walk in lattices, for example, in a model of turbu- 
lence which describes a transport in a Boltzmann lattice gas [28]. One can 
constrain the jumping size either by taking into account that the Brownian 
particle actually possesses a finite velocity, i.e. to substitute the Levy fiights 
by the Levy walks, or by introducing some cutoff in the jumping size prob- 
ability distribution. As a result, the variance becomes finite and, in the case 
of the mutually independent jumps, the random walk probability distribution 
converges to the Gaussian, according to the central limit theorem. The trun- 
cation of the Levy fiights can be accomplished either by a simple removing 
of the tail [29] or by multiplying the tail by some fast falling function. An 
obvious choice in this context is the exponential exp(— 7|x|) and the algebraic 
function where /3 > 2 — The former case was considered by Koponen 

in Ref. [30], where an analytic expression for the characteristic function was 
derived. The Levy fiights with exponential truncation serve as a model for phe- 
nomena in many fields, e.g. in turbulence [31], solar systems, economy. The 
distribution function of velocity and magnetic-field vector differences within 
solar wind can be reasonably fitted in this way [32] . In the framework of the 
economic research, the Levy process is a natural model of the financial assets 
fiow and the fractional equations are applied to characterize the dynamics 
of stock prices, where rare, non-Gaussian events are frequently encountered, 
in particular if the market exhibits high volatility. One can improve in this 
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a|x| 

Fig. 1. The distributions p{x,t) obtained from trajectory simulations for the pow- 
er-law cutoff (35) with fi' = 5.5 for t = 5, t = 10, t = 50 and t = 200 (from top to 
bottom). The distribution predicted by Eq.(30) is also presented and it coincides 

with the curve for t = 200. Other parameters are the following: 9 = 1, a = 0.1 and 
a = t-i/(2+9). 

way the Black-Scholes model, commonly used to price the options, which is 
restricted to the Gaussian distributions. In order to incorporate the finite- 
ness of the financial system to the fractional equations formalism by means 
of the exponential truncation of the Levy tails, models known as CGMY and 
KoBoL have been devised [33] . The non-Markovian fractional equation, which 
results from the CTRW model with the Levy distributed and exponentially 
truncated jump size, has been considered in Ref. [34]. The solutions do not 
exhibit a typical scaling at small time but they converge, asymptotically, to 
the stretched Gaussian which is predicted by the subdiffusive case of CTRW 
with the Gaussian step-size distribution. 

In Ref. [30], the following jumping size distribution, in the form of the Levy 
tail multiplied by the exponential factor, has been introduced: 

A(x) = Are-^l"l|xr'^-^ (32) 

and the characteristic function for the process has been evaluated. In the 
above formula 7 > and is the normalization constant. For the symmetric 
process, the normalized Fourier transform from A(x) is given by [30,34] 

A(A;) = -nT{fi) tan{7Tfi/2) [Y - {k^ + l^Y'^ cos(/i arctan(A;/7)] + 1. (33) 

TT 

We keep only the terms of the lowest order in \k\. The expansion of the ex- 
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Fig. 2. The time T needed to reach the distribution (30), as a function of 6, for the 
case of the algebraic cutoff (35) with the parameters n' = 5.5 and a = 0.1. 

pression (33) produces the following result 

X{k) ^ 1 + -/xr(/x)7''-Han(7r/x/2)(/x2 - /x)A;^ = 1 - K^k'^. (34) 



The diffusion process is then described by Eq.(29) with K'^ 



On the other hand, one can apply the power-law cutoff to the Levy tail. In 
Ref.[35] the fractional equation of the distributed order, with the constant 
diffusion coefficient, was introduced; it implies the cutoff in the form of the 
power-law function with the exponent 5 — /x. Since the equation involves both 
the fractional and the diffusion component, there is no simple scaling. In this 
paper, we assume the jumping size distribution in the form of the modified 
Levy tail: 



X{x) = < 



for \x\ < a 

ljL'a^'\x\~^'^^ for \x\ > (7, 



(35) 



where fj,' — fj, + P > 2 to ensure the existence of the second moment. The 
Fourier transform is given by [36] 

m = 1 - ^l^^C"^)" - - cos(7r/.72)] {\k\ar' + .... (36) 

In this case we have = /x'(T^/2(/x' — 2), provided we keep only the quadratic 
term. 
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100 200 300 400 500 

t 

Fig. 3. The shape parameter 6 of the distribution tail exp{—const\ax\^) as a function 
of time for the exponential cutoff (32). The curves correspond to the cases: 6 = 1, 

6 = and 6 = —0.5 (from top to the bottom). The other parameters: /x = 1.5 and 

7 = 1. 

In the following, we evaluate the probability density distribution p{x, t) for 
both forms of the Levy flight cutoff, exponential and algebraic, by means of the 
random walk trajectory simulations. The waiting time is sampled from Eq.(l) 
and the jumping size from either (32) or (35). Fig.l presents the power-law case 
for ^ = 1, which corresponds to the subdiffusion. The results are compared 
with the Kramers-Moyal limiting distribution (30), which is expected to be 
reached at large time. We observe a rapid convergence for small and medium 
|x|-values, whereas the tails reach the form (30) at about t = 200. Nevertheless, 
the shape of the tails is always stretched exponential ~ e^\){—const\ax\^) and 
the index 5 rises with time. The speed of convergence to the distribution (30) 
strongly depends on the parameter 9. In Fig. 2 we present the convergence 
time T as a function of 9. It is relatively short only for small |^|; for this 
case the kernel in the master equation (3) changes weakly with and the 
higher terms in the Kramers-Moyal expansion soon become negligible. T rises 
rapidly for the negative 9: the estimation presented in the figure suggests that 
the dependence T{9) is exponential, ~ exp(— 16^), which yields T ~ 10^ when 
9 approaches -1. The rapid growth of the time needed to reach convergence 
of the tails for the negative 9 may also be related to a specific shape of the 
distribution. The tails become flat for ^ < 0, and the asymptotics emerges 
first for very large On the other hand, the probability density to stay in 
the origin, p(0,t), is then infinite: we have p{x,t) ~ t^i'^+^)/i'^+^)\x\^ {\x\ <^ 1), 
according to Eq.(30). Note that in the case 6' = we obtain the usual result 
for diffusion: p{0,t) ~ l/Vi- 

Application of the exponential cutoff to the Levy tail produces similar prob- 
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0,6 0,8 1,0 1,2 1,4 

a|x| 

Fig. 4. The distributions p{x, t) at t = 100 and t = 200 obtained from trajectory 
simulations for the exponential cutoff with // = 1.5 and 7 = 2. These curves coincide 
and assume the numerically estimated shape ~ a\x\ exp(— 4.7(a|x|)^'^. The distribu- 
tion predicted by Eq.(30) is also presented. The parameter 6 = 1 and a = 

ability density distributions to those presented in Fig.l and they are also 
characterized by the stretched Gaussian tails. However, the convergence rate 
of the index 6 to the value 2 + 9, predicted by the solution (30), is smaller 
than for the power-law truncation because the kernel in Eq.(3) is steeper in 
this case. In Fig. 3 we present the dependence 6{t) for all kinds of the diffu- 
sion. Initially, the exponent rises fast but then it begins to stabilize and the 
curves approach the asymptotic values very slowly, especially for the superdif- 
fusive case of the negative 6. Discrepancies from Eq.(30) are more pronounced 
for 7 = 2. This case is presented in Fig. 4: though the rescaled distribution 
seems to be stabilized already for t = 100, its shape differs from (30) also for 
intermediate values of a|x|. 

The exponential asymptotics of the probability density distributions is guar- 
antied by the existence of the finite second moment. However, the time needed 
to reach that asymptotics can be so large [29] that it cannot be observed in any 
practical realizations of the process. In fact, the convergence rate to the normal 
distribution is governed by the third moment, according to the Berry-Esseen 
theorem [37] which refers to the sum of mutually independent variables, sam- 
pled from the same distribution. It states that if the third moment p is finite, 
then the deviation of a given distribution from the normal one is less than 
33 p/{4:a^^/n), where n the number of steps and a is the standard deviation. If 
p = oo, the convergence to the Gaussian may be problematic in practice. To 
demonstrate that the case of divergent third moment is exceptional also for 
our process, let us consider the power-law cutoff of the Levy tail, Eq.(35), with 
the parameter /x' = 2.5. The tail of the resulting random walk distribution, 
presented in Fig. 5, is no longer exponential but it assumes the power-law form 
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a|x| 

Fig. 5. The distributions p{x,t), plotted in the double logarithmic scale, obtained 
from trajectory simulations for the power-law cutoff with /x' = 2.5, 9 = 1 and 
a = 0.1. 

Ixl^*^ which persists to the largest times, numerically accessible. The param- 
eter 6 is constant and well determined for small time: 6 = 3.4. Moreover, it 
seems to be independent of 6. The presence of algebraic tail is not restricted 
to the case /i' < 3, for which the third moment is divergent; also for slightly 
larger values of fi' it is clearly visible. In the case /i' = 3.1 we find 5 = 4.0 at 
small times, whereas 6 = 4.5 for /j,' = 3.5. 

The probability distributions which possess algebraic tails with 3 < 5 < 4 are 
of interest in the economic research. It has been suggested that such power- 
laws in financial data arise when the trading behaviour is performed in an 
optimal way [38]. The stock market data seems to confirm that expectation. 
Extensive studies of the US indexes indicate the power-law form of the proba- 
bility distribution of stock price changes with 2.5 < S < 4 [39,40,41]. Moreover, 
the distributed order equation of Sokolov et al. [35] predicts the similar value: 
6 = 3.3. 



4 Conclusions 

We have solved the fractional equation which follows from the master equation 
for the jumping process in the diffusion approximation of small wave numbers. 
The jumping size distribution has the Levy form. The jumping rate depends on 
position and then the diffusion coefficient in the fractional equation is variable. 
We have considered the jumping rate in the algebraic form z/(x) = \x\~^, well 
suited for the diffusion on self-similar structures. The generalization to other 
dependences ul^x) is possible [7]. It has been demonstrated that the equation 
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is satisfied by the scahng formula which can be expressed in terms of the 
Fox function H22, when one neglects higher terms in the /c-expansion. In the 
limit /X — > 2, the solution reduces itself to the Fox function of the lower order 
and it exhibits the stretched-Gaussian asymptotic form (24). The exponent 
of that function, 1/a, is related to the higher moments and it cannot be 
uniquely determined in the diffusion approximation. The solution predicts all 
kinds of diffusion, both normal and anomalous, which are distinguished by the 
parameter 9. The diffusion equation can be solved exactly for the case fi = 2 
and that solution constitutes the Kramers- Moyal approximation of the master 
equation. The requirement that the fractional equation solution should agree 
with that result in the limits ji ^ 2 and t — > 00 allows us to find additional 
conditions for the Fox function coefficients. 

In the approximation of small wave numbers, the problem of truncated Levy 
flights coincides with that of the Gaussian jump sizes. We have applied two 
forms of the cutoff: the exponential and the algebraic ones to study the ran- 
dom walk process by the Monte Carlo method. In most cases, the probability 
density distributions converge with time to the Kramers-Moyal result which 
predicts a — 1/(2 + 6'). However, that convergence appears very slow if 9 is 
far away from 0, especially for ^ < 0. If the truncation function is steep, the 
distribution seems to assume the stabilized asymptotic shape which differs 
slightly from Eq.(30) (see Fig. 4) and this conclusion may indicate a limit of 
applicability of the Kramers-Moyal approximation. In other cases, form (30) 
is actually reached after a long time. For smaller time, the scaled distribution 
is time-dependent, in disagreement with the ansatz (7). However, the distribu- 
tion tail is always power-law and the exponent 5 depends on time very weakly 
(Fig. 3), compared to the function a{t). One can expect that in an experimen- 
tal situation, when the observation time is finite, the distributions are unable 
to converge to (30) and they may reveal the values of 6 smaller than 2-\- 9. 

The convergence of the distribution to (30) becomes problematic not only for 
very sharp cutoffs but, conversely, for the functions which possess a large, in 
particular infinite, third moment. Numerical simulations predict in this case 
the power-law asymptotics and no trace of the exponential tail could be found. 
The value of the exponent of power-law tail agrees with observations, e.g. for 
the financial data. 
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